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ABSTRACT 


A  quantitative  evaluation  of  the  limitations  of  the  radiation  boundary 
elements  in  the  finite  element  code  ATILA  [Ref.  1]  has  been  performed.  Five 
three-dimensional  models  were  employed,  each  representing  a  rigid  spherical 
solid  surrounded  by  water.  Monopolar,  dipolar  and  quadrupolar  incident 
spherical  waves  were  introduced  and  the  corresponding  scattered  waves  were 
computed  using  the  ATILA  code  and  an  exact  analytical  solution. 

The  dimensionless  parameters  that  characterize  the  problem  are  ka,  kL, 
and  kR  where  k  is  the  wavenumber  of  sound  in  water,  a  is  the  radius  of  the 
scatterer,  R  is  the  outer  fluid  mesh  radius,  and  L  is  the  thickness  of  the  fluid 
layer.  The  range  of  values  investigated  were  kR=1.5,  2.5,  4.0,  ka=0.5,  1.0,  2.0 
and  kL=0.5,  1.0. 

For  axially  symmetric  incident  fields,  the  maximum  normalized  errors 
occurred  at  the  poles  and  were  9%,  12%,  and  6%,  respectively.  Furthermore, 
the  errors  for  monopolar  and  dipolar  incident  fields  were  strongly  influenced  by 
the  location  of  the  radiation  boundary  (kR),  less  so  by  the  scatterer’s  radius  (ka); 
specifically,  the  error  decreases  with  increasing  kR  and/or  ka.  The  errors  for 
quadrupolar  incident  fields  do  not  exhibit  any  significant  dependence  on  kR  or 
ka.  The  errors  for  all  the  axially  symmetric  incident  fields  were  not  affected  by 
variations  of  the  element’s  size  (kL).  For  non-axially  symmetric  incident  fields, 
the  maximum  deviation  occurred  at  the  equatorial  points  and  was  less  than 
5.5%. 

Further  investigation  using  a  two-dimensional  model  is  proposed  in  order 
to  determine  the  range  of  values  of  ka,  kL,  and  kR  which  will  result  in  negligibly 
small  errors. 
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I.  INTRODUCTION 


The  investigation  described  in  this  thesis  is  part  of  an  ongoing  research 
project  in  the  numerical  modeling  of  arbitrary  densely  packed,  random  volumetric 
active  sonar  arrays  [Ref.  1].  The  technique  employed  is  an  extension  of  the  so- 
called  T-matrix  method,  which  has  been  applied  to  other  scattering  problems 
[Ref.  2].  This  method  rigorously  accounts  for  multiple  scattering  to  all  orders.  To 
apply  the  T-matrix  method  to  the  problem  of  active  sonar  array  performance 
prediction,  it  is  necessary  to  compute  the  radiating  and  scattering  properties  of  a 
single  array  element  in  a  free  field  environment.  To  accomplish  this  for  a  real 
transducer,  the  ATILA  [Ref.  3]  finite  element  code  is  employed. 

ATILA  is  a  French  finite  element  code  especially  developed  for  the 
analysis  of  underwater  acoustic  transducers.  This  code  was  written  by  engineers 
at  the  Institut  Superieur  d'Electronique  du  Nord  (ISEN),  Lille,  France  and  is  in 
use  by  U.S.  scientists  within  and  outside  the  Navy  working  on  U.S.  Navy- 
sponsored  research. 

Computation  of  the  radiating  and  scattering  properties  of  a  sonar 
transducer  using  ATILA  involves  building  a  finite-element  model  representing  the 
transducer  and  surrounding  it  by  a  finite-element  mesh  representing  water,  which 
is  terminated  by  so-called  “radiation  boundary  elements”.  The  radiation 
boundary  elements  are  intended  to  “absorb”  all  incident  acoustic  radiation.  In 
practice,  they  perform  this  function  less  than  perfectly.  This  has  a  direct  effect 
on  the  computation  of  a  transducer’s  radiating  and  scattering  properties. 

The  present  investigation  focuses  on  the  evaluation  of  ATILA's  radiation 
boundary  elements.  For  this  purpose,  a  three  dimensional  spherical  fluid  mesh, 
surrounding  a  spherical  rigid  body  and  terminated  by  radiation  boundary 
elements  is  studied.  An  incident  spherical  wave  consisting  of  a  single 
component  of  a  multipolar  acoustic  field  strikes  the  sphere,  and  the  resulting 
scattered  wave  field  is  calculated  using  ATILA.  The  amounts  of  acoustic  field  in 
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the  intended  outgoing  component  and  in  other  multipolar  components  are  then 
analyzed.  The  influence  on  the  results  of  the  fluid  mesh  element  size,  fluid  mesh 
inner  radius,  and  fluid  mesh  outer  radius  relative  to  the  acoustic  wavelength  are 
investigated  and  quantitative  guidelines  developed  in  order  to  minimize  the  effect 
of  imperfect  absorption. 

The  remainder  of  this  thesis  is  divided  into  four  chapters.  Chapter  II 
describes  the  theory  involved  in  the  finite  element  analysis,  the  possible  types  of 
analyses  that  can  be  performed,  particularly  the  harmonic  analysis  of  a  radiating 
spherical  structure  excited  by  an  incident  spherical  wave. 

Chapter  III  describes  the  three  dimensional  spherical  models  which  were 
employed.  Chapter  IV  presents  and  discusses  the  scattering  results  for  different 
incident  multipolar  components  and  the  influence  of  element  size  and  mesh 
inner  and  outer  radius.  Chapter  V  concludes  the  thesis.  Appendix  A  contains 
the  C  program  used  to  analytically  calculate  the  pressure  scattered  from  the 
spherical  boundary  for  pressure  release  and  rigid  surface  boundary  conditions. 
Appendix  B  presents  the  FORTRAN  code  used  in  ATILA  for  the  generation  of 
the  spherical  incident  pressure  wave. 
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11.  THEORY 


A.  GENERAL  PRESENTATION  OF  ATILA  VERSION  5.03 

The  finite  element  method  is  a  technique  that  provides  numerical 
solutions  for  boundary  value  problems  and  field  calculations.  [Refs.  4, 5, 6, 7, 8] 
ATILA  is  a  finite  element  code  developed  specifically  for  the  analysis  of  sonar 
transducers. 

The  ATILA  code  is  able  to  perform: 

1 .  elastic  or  piezoelectric  structures  modeling, 

2.  magnetostrictive  structures  modeling, 

3.  periodic  structures  modeling, 

following  a  general  formulation  shown  in  Figure  1  (after  [Ref.  3]),  to  model  an 
undenwater  acoustic  transducer. 


Figure  1.  ATILA  general  formulation.  After  Ref.  [3] 
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The  ATII_A  code  is  based  on  the  separation  of  the  physical  problem  under 
consideration  into  a  discrete  number  of  elements  which  are  in  turn  described  by 
their  nodes,  in  a  given  order.  For  each  node  there  are  a  number  of  degrees  of 
freedom  (d.o.f.)  that  can  be  specified  using  certain  boundary  conditions.  The 
elements,  nodes,  and  the  specific  node-numbering  order,  is  referred  to  as  the 
topology  of  the  problem. 

An  ATILA  job  organization  is  carried  out  in  several  steps,  as  follows: 

1 .  Model  definition. 

2.  Mesh  generation. 

3.  Data  file  preparation. 

4.  Running  a  job. 

5.  Result  file  postprocessing. 

First  of  all,  the  type  of  analysis  to  be  performed  has  to  be  specified,  for 
example,  harmonic  analysis  of  a  solid  structure  excited  by  an  impinging  wave. 
The  type  of  elements  that  are  required  to  describe  the  fluid  and  structure  domain 
are  then  chosen.  ATILA  includes  a  library  of  46  different  elements  to  model 
composite,  piezoelectric,  fluid,  magnetostrictive,  coupling  FEM-BEM,  interface, 
and  radiation  dampers. 

The  mesh  generation  procedure  includes  the  assignment  of  coordinates 
for  each  node  and  the  node-numbering  order  for  each  element.  Throughout  this 
procedure  the  whole  physical  problem  is  discretized  into  elements  which  allow 
the  representation  and  modeling  of  many  different  geometrical  shapes  and  lines, 
as,  for  example,  PRIS15F,  a  fifteen-node  isoparametric  triangular  base  prism 
used  to  model  homogenous  fluid  media,  or  TRIA06R,  a  six  node  isoparametric 
triangular  element  to  prescribe  a  monopolar,  dipolar,  or  multipolar  radiation 
condition. 

One  of  the  facilities  of  the  ATILA  code  is  the  MOSAIQUE  preprocessor. 
This  routine  enables  the  use  of  super-elements  and  generates  all  the  necessary 
elements  and  node  data  for  ATI  A. 


4 


The  data  file  includes  all  the  neccessary  input  data  such  as  the  type  of 
analysis,  the  material  properties,  the  node  coordinates,  the  elements,  the 
boundary  conditions,  the  loading  data,  and  possibly  the  plane  wave  data,  to  carry 
out  the  analysis. 

Running  an  ATILA  job  involves  calling  up  the  subroutines  to  compute 
elementary  matrices,  solve  equations  and  display  the  results. 

The  available  results  file,  if  desired,  can  be  postprocessed  to  create 
graphic  displays  of  the  structure,  contours  of  constant  value  for  potentials, 
pressures,  and  displacements.  A  simple  flow  chart  of  an  ATILA  job  is  shown  in 
Figure  2  (after  [Ref.  3]). 


Figure  2.  Flow  chart  of  an  ATILA  job.  After  Ref.[3] 
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B.  ATILA  FINITE  ELEMENT  COMPUTATION  FOR  RADIATING  AND 
SCATTERING  PROBLEMS 

1.  General  Elastic  or  Piezoelectric  Structures  Modeling 
A  large  number  of  analyses  can  be  performed.  These  are  based  upon  the 
complete  set  of  equations  of  elasticity  in  the  structure,  the  Helmholtz  equation  in 
the  fluid,  and  Poisson's  equation  in  the  elastic  or  piezoelectric  material. 
Appropriate  radiation  boundary  conditions  are  applied  on  the  surface  which 
surrounds  the  fluid  domain. 

The  unknown  quantities  are  the  nodal  values  of  the  displacement  field  U 
in  the  whole  structure,  the  electric  potential  O  in  the  piezoelectric  material,  and 
the  pressure  P_in  the  fluid.  The  system  of  equations  is  written  in  matrix  form: 


■([K,J-(o2[M])  [K,^]  -[L] 

y‘ 

F 

[K^l  [0] 

= 

"9 

-pVco2[Lr  [Of  ([H]-a)2[Mi])_ 

p 

pc^  V 

where: 

U:  vector  of  the  nodal  values  of  the  components  of  the  displacement  field. 
O:  vector  of  the  nodal  values  of  the  electric  potential. 

P:  vector  of  the  nodal  values  of  the  pressure  field. 

F:  vector  of  the  nodal  values  of  the  applied  forces, 
g:  vector  of  the  nodal  values  of  the  electric  charges. 

vector  of  the  nodal  values  of  the  integrated  normal  derivative  of  the 
pressure  on  the  surface  boundary  S. 

[Kuj:  stiffness  matrix. 

[KuJ :  piezoelectric  matrix. 

[K,j^]:  dielectric  matrix. 

[M]:  consistent  mass  matrix. 

[H]:  fluid  (pseudo-)stiffness  matrix. 


6 


[Mi]:  consistent  (pseudo-)fluid  mass  matrix. 

[L]:  coupling  matrix  at  the  fluid  structure  interface  (connectivity  matrix). 

[0]:  zero  matrix. 

co;  angular  frequency. 

p:  fluid  density. 

c:  fluid  sound  speed. 

T:  means  transposed. 

ATILA  is  able  to  perform: 

1.  Static  analysis  of  an  elastic,  piezoelectric,  or  hydroelastic  system, 
where  the  displacement  field  and/or  the  electric  potential  or  the  pressure  fields 
are  required. 

2.  Modal  analysis  of  an  elastic,  piezoelectric,  hydroelastic  system, 
where  the  eigenfrequencies,  the  resonance  and  antiresonance  frequencies,  and 
the  normal  modes  are  computed. 

3.  Harmonic  analysis  of  a  driven  elastic  or  piezoelectric  structure,  or 
the  scattering  of  an  arbitrary  incident  wave  by  an  elastic  or  piezoelectric 
structure. 

A  scattered  wave  analysis  of  a  spherical  pressure  wave  incident  on  a  solid 
structure  is  required  to  investigate  the  performance  of  the  radiation  boundary 
elements,  and  is  presented  in  the  following  section. 

2.  Harmonic  Analysis  of  a  Solid  Structure  Excited  by  an 
Impinging  Wave 

In  this  type  of  analysis,  the  real  and  imaginary  parts  of  the  pressure  field 
(P),  the  displacement  field  (U),  the  potential  (^,  and  the  electric  current  (icoQ) 
are  computed.  The  pressure  (P)  and  the  flux  ,  are  split  into  incident  and 
scattered  parts.  The  normal  derivative  of  the  incident  pressure  on  the  surface 

boundary  S  is  written  'Fj=[D] — where  [D]  is  a  matrix  obtained  by  assembling  the 

damping  elements  on  the  surface,  provided  by  the  code. 
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The  incident  pressure  field  is  provided  via  a  FORTRAN  function 

"INCPRE  (x,y,z,k)",  included  at  the  end  of  the  main  program  by  the  user  as 
shown  in  Appendix  B,  after  [Ref  2],  Utilization  of  a  user-provided  incident 
pressure  allows  the  excitation  of  the  structure  with  an  incident  wave  of  any  kind, 
while  the  default  function  provided  with  ATI  LA  generates  a  plane  wave  e*'" 
traveling  from  the  positive  to  the  negative  direction  along  the  Ox  axis  (ej“*  time 
harmonic  dependence  is  assumed). 

By  assigning  a  specific  entry  in  the  ATI  LA  code,  we  are  able  to  compute 
the  total  pressure  or  just  the  scattered  one,  utilizing  the  "TOTAL"  or 
"SCATTERED"  attribute,  respectively. 

The  equations  of  elasticity  in  the  structure,  the  Helmoltz  equation  in  the 
fluid,  and  Poisson's  equation  in  the  solid  material  are  written  in  matrix  form  : 


{[K,J-a)2[M]) 

[Ku*r 

-[Lf 


where: 

Pes :  vector  of  the  nodal  values  of  the  elastic  scattered  pressure  field, 

Pi :  vector  of  the  nodal  values  of  the  incident  pressure  field, 

[G]:  frequency  dependent  complex  linear  matrix, 

[D]:  matrix  of  the  damping  elements. 

Furthermore,  the  internal  losses  of  the  material  used  can  be  taken  into 
account  throughout  a  specified  program  entry  "SKYLINE  COMPLEX"  and  when 
no  losses  are  required,  with  "SKYLINE  REAL",  respectively. 


[Ku*] 

[0]^ 


[L] 

[0] 

[H] 


VpVco^ 


IMd 


y 

d) 


2  2  2 
p^C^O)^ 


[G]Pe 


F-[L]Pi 

rQ 
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111.  THREE  DIMENSIONAL  SPHERICAL  MODELS 


A.  INTRODUCTION 

In  order  to  evaluate  ATILA’s  radiation  boundary  elements,  we  first 
developed  a  family  of  three-dimensional  spherical  models.  These  models 
simulate  a  spherical  rigid  structure  surrounded  by  an  infinite  fluid  environment. 

A  total  of  five  models  were  employed.  All  models  are  composed  of  an 
inner  spherical  boundary  representing  a  rigid  solid,  several  concentric  spherical 
fluids  layers,  and  an  outer  spherical  boundary  composed  of  radiation  boundary 
elements  (dampers),  each  representing  a  semi-infinite  fluid  region. 

B.  DESCRIPTION  OF  MODELS 

The  models  developed  are  distinguished  by  their  values  of  the  inner 
radius  a,  the  number  and  thickness  L  of  each  fluid  layer,  and  the  outer  radius  R. 
The  appropriate  dimensionless  forms  for  these  parameters  are  ka,  kL,  kR,  where 

k=—  is  the  acoustic  wavenumber.  Table  1  lists  the  properties  of  each  model 
c 

used  in  the  investigation. 

Model  1  had  already  been  used  in  the  calculation  of  the  transition  matrix 
for  the  scattering  of  acoustic  waves  from  a  thin  elastic  spherical  shell  [Ref.  2]. 
This  model  contains  four  fluid  layers  of  two  different  thicknesses.  To  separately 
investigate  the  influence  of  each  of  the  dimensionless  parameters  (ka,  kR,  kL), 
four  additional  models  were  created. 

Model  2  serves  as  a  reference  model.  Each  of  models  3  through  5  is 
distinguished  from  model  2  in  that  the  value  of  only  one  of  ka,  kL,  kR  is  changed. 


9 


Model  1 

Model  2 

Model  3 

Model  4 

Model  5 

Original 

Basic 

ka 

kL 

kR 

Wavenumber  k  in  meters  ’ 

2 

1 

1 

1 

1 

Radius  a  in  meters 

0.5 

0.5 

1.0 

0.5 

0.5 

Scatterers  Radius  ka 

1.0 

0.5 

1.0 

0.5 

0.5 

Number  of  Fluid  Layers 

2  of 

4  of 

3  of 

2  of 

2  of 

L  =  0.25 

L  =  0.5 

L  =  0.5 

L=  1.0 

L  =  0.5 

2  of 

1“ 

II 

O 

cn 

Element’s  Size  kL 

1.0 

0.5 

0.5 

1.0 

0.5 

Radiation  Boundary  kR 

4.0 

2.5 

2.5 

2.5 

1-5 

Table  1 .  Properties  of  each  model  used  in  the  investigation. 

The  solid  structure  is  modeled  by  specifying  a  zero-flux  or  boundary 
condition  on  the  scatterer  surface  (interface  between  solid  and  fluid  domain). 
This  is  the  default  condition. 

For  radiation  problems,  the  fluid  mesh  outer  limit  must  be  spherical. 
Therefore,  the  surrounding  solid  structure  fluid  is  modeled  with  a  spherical 
surface  of  dimensionless  radius  kR  =  2.5  or  1.5,  for  models  2  through  5. 

The  fluid  (water)  is  simulated  by  assuming  the  following  properties: 

1.  E=0.222*l0’°Pa  (Young's  modulus) 

2.  p=0.1*10‘*kg/m®  (Density) 

3.  v=0.0  (Poisson's  Ratio) 

The  whole  solid  structure  and  fluid  mesh  was  constructed  using  fifteen- 
node,  three-dimensional  isoparametric  right  triangular  prismatic  elements 
(PRIS15)  from  the  ATILA  finite  element  library  [Ref.3]. 

ATILA  offers  monopolar,  dipolar,  and  multipolar  radiation  boundary 
elements.  Multipolar  damping  elements  were  selected  to  terminate  the  fluid 
mesh  outer  radius  and  the  appropriate  condition  was  set  in  the  data  input  file. 
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The  multipolar  damping  element  used  is  a  six-node  isoparametric  triangular 
element  (TRIA06R)  and  it  has  to  be  spherical  and  attached  to  the  outer  surface 
of  the  fluid  mesh. 

The  topology  of  the  prismatic  (PRIS15)  and  of  the  triangular  (TRIA06R) 
elements  are  shown  in  Figures  3  and  4,  respectively.  The  numbers  represent 
the  nodes  and  the  order  of  numbering  required  for  each  element. 

C.  MESH  GENERATION 

The  major  limitation  of  the  three-dimensional  finite  element  model  is  the 
number  of  degrees  of  freedom  (D.O.F.)  available  (displacement  along  the 
coordinate  axes,  rotation  along  the  coordinate  axes,  pressure,  electric  or 
magnetic  potential,  and  magnetic  excitation  currents).  For  this  reason  the 
models  developed  were  limited  by  the  number  of  nodes  and  elements  that  were 
allowed. 

Specifically,  the  number  of  nodes  employed  in  the  models  ranged  from 
1546  to  2346  for  the  structure  and  fluid  mesh,  and  the  total  number  of  elerhents 
ranged  from  216  to  360.  Of  these,  144  to  288  are  fluid  elements  and  72  are 
damping  elements.  The  fluid  mesh  is  arranged  in  successive  layers  of  various 
thicknesses.  The  original  fluid  mesh  was  divided  into  four  layers,  two  of  them 
with  0.25m  thickness  and  the  other  two  fluid  layers  with  0.5m  thickness. 

Figure  5  shows  a  Mercator  projection  of  the  nodes  and  elements  on  the 
inner  surface  of  the  fluid.  The  node  numbers  and  selected  element  numbers 
corresponding  to  this  layer  are  given. 

For  the  mesh  generation  procedure  all  the  nodes  were  specified  using 
spherical  coordinates  with  the  origin  at  the  center  of  the  spherical  structure.  The 
mesh  spacing  was  always  less  than  a  quarter  of  the  acoustic  wavelength  used. 

The  aspect  ratio  of  each  element  according  to  the  reference  manual 
should  be  less  than  or  equal  to  3;  in  our  models  it  is  1  to  4. 
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Figure  4.  TRIA06R  Element  Topology.  After  Ref. [3] 
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Figure  5.  Inner  Surface  Fluid  Layer  Model,  Mercator  Projection. 
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The  internal  angles  of  our  elements  were  modeled  between  30  and  90 
degrees  while  the  reference  manual  states  that  these  should  lie  between  30  and 
100  degrees. 

The  mesh  was  built  in  such  a  way  that  adjacent  elements,  like  elements  1 
and  2  shown  in  Figure  5,  were  able  to  share  common  nodes.  The  top  side 
exploded  mesh  view  of  the  three  dimensional  model  created  by  ATI  LA  DEPL 
program  is  shown  in  Figure  6. 

The  types  of  isoparametric  elements  used  in  the  mesh  generation  are 
described  in  [Ref.  3]  and  listed  in  the  following  Table  2. 


Region 

Element 

Geometry 

Fluid 

PRIS15F 

15-node  triangular 

prism 

Radiation  Damping 

TRIA06R 

6-node  triangular 

Table  2.  Isoparametric  Elements  Used  in  The  Three  Dimensional  Mesh 
Generation 
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Figure  6.  Three  Dimensional  Spherical  Model.  Topside  View 
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IV.  RESULTS 


A.  EXACT  ANALYTICAL  RESULTS  FOR  SCATTERED  PRESSURE  WAVE 

The  linearized  homogenous  wave  equation  for  the  propagation  of  sound 
in  ideal  (nonviscous)  fluids  and  the  time-independent,  lossless  Helmholtz 
equation  for  a  pressure  wave  with  time-harmonic  dependence  p(t,r)=p(r)*e''“* ,  at  a 
position  r=(r,0,(p)  and  time  (t),  is  given  by  [Refs.  9,  10] ; 

V^p(t,[)-4-  =  0  =>  v2p(r)  +  k2p{r)  =  0  (1 ) 

The  solution  of  the  Helmoltz  Equation  (1)  for  the  incident  (Pm)  and 
scattered  (Pjc)  pressure  fields  in  the  spherical  coordinate  system  shown  in 
Figure  7  is  given  by  Equations  2  and  3: 


Figure  7.  The  Spherical  Coordinates  (r.G.cp). 


(2) 


P|n(r.e.(p)=  I  Z  Pn^1-h};')(kr)-Y;;'(e,<t)) 

n  =  Om  =  -n 

00  n  /o\ 

Pso(r.e.v)=  Z  Z  Pnm2''n'(''^)C(®'« 

n  =  Om  =  -n 


where: 

Pnmi  and  Pnm2  3''e  the  amplitudes  of  the  n,m  components  of  the  incident 
and  scattered  waves,  respectively, 

h<,^^(kr)  =  the  nth  order  spherical  Hankel  function  of  the  first  kind, 
h[,^>(kr)  =  the  nth  order  spherical  Hankel  function  of  the  second  kind, 

=  the  spherical  harmonic  of  order  n,m,  which  is  related  to  the 
associated  Legendre  polynomial  by  the  equation 

Application  of  boundary  conditions  on  the  inner  surface  of  the  spherical 
fluid  mesh  for  vanishing  of  the  pressure  or  its  normal  derivative,  provides 
solutions  for  the  scattered  wave  for  the  case  of  a  pressure  release  or  rigid 
surface. 

1.  Analytical  Expression  for  the  Scattered  Pressure  from  a 
Pressure  Release  Spherical  Surface 

The  following  equations  apply  for  the  n,m  component: 


Pin(r.0.9)=Y^'(e,(p)-hn(^\kr) 


(4) 


Psc(r.0.<P)=[-h[,''^ka)  /  h(,2)(ka)]  •  Y;j'(e,(p)  •  h(,2)(kr)  (5) 
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where: 


a  =  the  radius  of  the  spherical  shell, 
k  =  the  wavenumber  of  sound, 

2.  Analytical  Expression  for  the  Scattered  Pressure  from  a  Rigid 
Spherical  Surface 


The  following  equations  apply; 


Pjn(r,0.(p)=Yj^{0.cp)-hJ,'')(kr),  (6) 

Psc('-.e.9)=[-(hSi'’^(ka))  /(h}f)(ka))]  Yj^(0.(p)  hJ,2)(,,r)_  (7) 

where  (h<,’^(ka))' and  (hjf>(ka))'are  the  spatial  derivatives  of  hf,^>(kr)  and 
evaluated  at  the  spherical  surface. 

Appendix  A  contains  a  C  program  used  to  calculate  analytically  the  values 
of  the  pressure  scattered/radiated  from  a  vibrating,  spherically  symmetric 
surface,  for  the  above  boundary  conditions  [Ref.  11]. 

B.  NORMALIZED  SCATTERED  PRESSURE  DEVIATION  COMPUTED  BY 

ATILA 

Using  the  requried  ATILA  input  data  file,  the  pressure  scattered  by  a  rigid 
spherical  surface  due  to  an  incident  spherical  harmonic  wave  was  computed. 
The  results  were  compared  with  an  exact  analytical  solution  computed  using  the 
C  program  described  in  Appendix  A.  Accordingly,  the  limitations  of  the  boundary 
elements  were  obtained. 

Since  the  scattered  pressure  is  proportional  to  the  incident  pressure  and 
has  angular  dependence,  the  error  In  the  ATILA  results  at  each  finite-element 
node  was  quantified  by  normalizing  the  deviation  from  the  exact  value  by  the 
maximum  magnitude  of  the  scattered  pressure  at  the  same  radius.  This 
represents  the  normalized  scattered  pressure  deviation  computed  by  ATILA. 
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C.  RESULTS  FOR  THE  THREE  DIMENSIONAL  SPHERICAL  MODEL 

(MODEL  1) 

The  following  Figures  8  and  9  depict  the  results  of  the  calculation  for  the 
original  three-dimensional  spherical  model  (model  1)  that  had  been  used  in  the 
previous  investigation  [Ref.  2],  In  this  model,  the  fluid  is  divided  into  four  layers 
of  thickness: 

1 .  Layer  a:  0.25  m 

2.  Layer  b:  0.25  m 

3.  Layer  c:  0.5  m 

4.  Layer  d:  0.5  m 

The  radius  of  the  scattering  (inner)  surface  (a)  is  0.5m:  the  radius  of  the  radiation 
boundary  is  2.0m.  The  frequency  (f)  is  474Hz;  the  corresponding  wavenumber 

of  sound  in  the  water  is  k=— =2m"l 

c 

Figure  8,  from  top  to  bottom,  presents  the  normalized  scattered  pressure 
deviation  computed  by  ATILA  using  multipolar  dampers  for  monopolar  and 
dipolar  incident  pressure  fields,  versus  the  dimensionless  radius  (kr)  from  the 
center  of  the  structure.  Figure  9  presents  the  normalized  scattered  pressure 
deviation  for  the  quadrupolar  incident  pressure  field. 

The  following  Table  3  provides  the  normalized  maximum  error  in 
percentages  in  the  middle  of  each  fluid  layer  and  the  angular  location  of  that 
error  for  the  monopolar,  dipolar,  and  quadrupolar  incident  fields. 

Nodes  located  at  the  poles  are  points  with  coordinate  spherical  angles:  0 
(polar)  =  000,180  and  (p  (azimuthal)  =  000,  degrees.  Nodes  located  on  the 
equator  are  points  with  coordinate  spherical  angles:  0  (polar)  =  090  and  9 
(azimuthal)  =  000,090,180,270  degrees.  Near  the  equator,  nodes  are  points  with 
coordinate  spherical  angles:  0  (polar)  =  075,105  and  9  (azimuthal)  = 
090,180,270  degrees. 
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NORMALIZED  SCATTERED  PRESSURE  DEVIATION 
COMPUTED  BY  ATILA  USING  MULTIPOLAR  DAMPERS 


MONOPOLAR  RADIATION 
n=0  /  m=0 


DIPOLAR  RADIATION 
n=1  /  m=0 


DIPOLAR  RADIATION 
n=1  /  m=±1 


Figure  8.  Normalized  scattered  pressure  deviation  computed  by  ATILA  versus  kr 
for  monopolar  and  dipolar  incident  fields. 
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NORMALIZED  SCATTERED  PRESSURE  DEVIATION  COMPUTED  BY 
ATILA  USING  MULTIPOLAR  DAMPERS 


QUADRUPOLAR 
RADIATION 
n=2  /  m=0 


0.06 
0.05 
0.04 
0.03 
0.02 
0.01 
0 

0  1  2  3  4  5 


QUADRUPOLAR 
RADIATION 
n=2  /  m=±1 


0.06 
0.05 
0.04 
0.03 
0.02 
0.01 
0 

0  1  2  3  4  5 

kr 


QUADRUPOLAR 
RADIATION 
n=2  /  m=±2 


Figure  9.  Normalized  scattered  pressure  deviation  computed  by  ATILA  versus  kr 
for  quadrupolar  incident  field. 
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Type  of 

Field 

Scatterer 

Surface 

ka=1.0 

Center 

of 

Layer  a 

kr=1.25 

Center  of 

Layer  b 

kr=1.75 

Center  of 

Layer  c 

kr=2.5 

Center  of 

Layer  d 

kr=3.5 

Maximum 

Normalized 

Error 

kr=1.5 

Monopolar 

2.6 

3.8 

5.7 

7.8 

8.6  kr=3.0 

4.6 

m=n=0 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Dipolar 

2.2 

3.0 

9.2 

8.5 

11.8kr=3.0 

3.9 

n=1,  m=0 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Dipolar 

1.8 

2.2 

3.6 

5.1 

4.0 

5.7  kr=3.0 

3.0 

n=1,  m=±1 

Equator 

Equator 

Equator 

Equator 

Equator 

Equator 

Equator 

Quadrupolar 

0.3 

0.4 

0.8 

2.7 

5.3 

5.8  kr=3.0 

0.6 

n=2,  m=0 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Quadrupolar 

0.3 

0.4 

0.5 

0.9 

3.7 

3.6  kr=3.5 

0.5 

n=2,  m=±1 

Equator 

Equator 

Equator 

Equator 

Equator 

Equator 

Equator 

Quadrupolar 

0.3 

0.5 

0.9 

1.9 

3.6 

3.6  kr=3.5 

0.5 

n=2,  m=+2 

Near 

Near 

Near 

Near 

Near 

Near 

Near 

Equator 

Equator 

Equator 

Equator 

Equator 

Equator 

Equator 

Table  3.  Percent  normalized  scattered  pressure  deviation  and  location  of  the 
maximum  error  for  monopolar,  dipolar  and  quadrupolar  incident  fields. 

From  the  above  table  and  Figures  8  and  9,  it  can  be  concluded  that: 

1.  The  greatest  normalized  errors  appear  for  the  following  type  of 

fields: 

a.  Monopolar,  n=m=0:  8.6%  at  kr=3.0 

b.  Dipolar,  n=1 ,  m=0:  1 1 .8%  at  kr=3.0 

c.  Quadrupolar,  n=2,  m=0:  5.8%  at  kr=3.0 

2.  For  the  above  axisymetrical  incident  pressure  fields  the 
corresponding  location  of  the  maximum  error  points  is  close  to  the  poles  and 
exactly  on  the  poles  for  the  scatterer’s  surface  and  on  the  poles  for  each  layer, 
respectively. 
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3.  For  the  cases,  where  the  maximum  error  points  appear  on  the 
equator  or  near  the  equator  (dipolar  n=1,  m=±1,  quadrupolar  n=2,  m=±1),  the 
location  of  the  minimum  error  points  is  on  the  poles. 

4.  If  the  maximum  error  points  are  ignored  for  the  monopolar  and 
dipolar  fields  then  the  maximum  normalized  error  is  always  less  than  4%. 

5.  If  the  maximum  error  points  are  ignored  for  the  quadrupolar  field 
then  the  maximum  normalized  error  is  always  less  than  3.5%. 

6.  The  normalized  scattered  pressure  deviation  computed  by  ATILA 
as  it  is  presented  in  Figures  8  and  9  appears  to  increase  when  moving  further 
away  from  the  acoustic  center. 

The  above  results  for  the  monopolar,  dipolar  and  quadrupolar  incident 
fields  indicate  that,  for  a  given  value  of  n,  the  maximum  normalized  errors  occur 
for  the  axially  symmetric  type  of  fields  (i.e.,  for  m  =  0).  Hence,  in  investigating 
the  influence  of  the  fluid  mesh  inner  radius  (ka),  the  element’s  size  (kL),  and  the 
fluid  mesh  outer  radius  (kR)  on  the  results,  only  the  axially  symmetric  incident 
waves  were  analyzed. 

D.  INFLUENCE  OF  FLUID  MESH  INNER  RADIUS 

In  order  to  evaluate  the  influence  of  fluid  mesh  inner  radius  on  the  results, 
models  2  and  3  were  used.  Recall  that  model  2  is  divided  into  four  equal  fluid 
layers  of  thickness  L=0.5m  and  the  scatterer’s  radius  is  0.5m.  Model  3  is  divided 
into  three  equal  fluid  layers  of  thickness  L=0.5m  and  scatterer’s  radius  1  .Om.  For 
both  cases,  the  radiation  boundary  and  element’s  size  satisfy  kR=2.5  and 
kL=0.5,  with  k  =  1.0m'^  and  R=2.5m. 

Figure  10  presents  the  influence  of  the  scatterer’s  radius  (ka)  on  the 
normalized  scattered  pressure  deviation,  versus  kr.  The  monopolar  (a,b),  dipolar 
(c,d)  and  quadrupolar  (e,f)  fields  are  shown  from  top  to  bottom.  On  the  left  side 
the  scatterer  radius  is  ka=0.5  and  on  the  right  side  the  scatterer  radius  is  ka=1 .0. 
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4-EQUAL  FLUID  LAYERS  0.5M 


3-EQUAL  FLUID  LAYERS  0.5M 


Table  4  summarizes  the  results  for  the  maximum  error  in  percent  in  the 
middle  of  each  layer  and  the  node  location  for  the  monopolar,  dipolar  and 
quadrupolar  incident  fields  when  the  scatterer  radius  ka  varies  and  kL  and  kR  do 
not. 


Type  of 

Field 

Scatterer 

Surface 

Center  of 

Layer  a 

kr=0.75 

Center  of 

Layer  b 

kr=1.25 

Center  of 

Layer  c 

kr=1.75 

Center  of 

Layer  d 

kr=2.25 

Maximum 

Normalized 

Error 

kr=1.5 

Monopolar 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

n=m=0 

ka=1.0 

2.4 

— 

3.5 

5.8 

5.7 

7.0  kr=2.0 

4.6 

ka=0.5 

1.7 

2.9 

4.9 

6.9 

5.9 

7.8  kr=2.0 

5.7 

Dipolar 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

n=1,  m=0 

ka=1.0 

1.5 

— 

1.9 

3.3 

3.8 

4.2  kr=2.0 

2.5 

ka=0.5 

0.5 

0.6 

1.8 

3.0 

3.9 

4.0  kr=2.0 

2.2 

Quadrupolar 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

n=2,  m=0 

ka=1.0 

0.4 

0.6 

1.4 

2.8 

4.1  kr=2.5 

0.9 

ka=0.5 

0.2 

0.3 

0.5 

1.2 

2.6 

4.1  kr=2.5 

0.7 

Table  4.  Percent  normalized  scattered  pressure  deviation  and  location  of  the 
maximum  error  for  monopolar,  dipolar  and  quadrupolar  incident  pressure  field 
when  ka=1.0  and  ka=0.5,  radiation  boundary  kR=2.5,  element’s  size  kL=0.5. 

From  this  table  and  Figure  10  it  is  observed  that: 

1 .  For  the  monopolar  field: 

a.  The  normalized  error  at  a  given  value  of  kR  decreases  as  ka 
is  increased  from  0.5  to  1 .0. 

b.  The  maximum  error  decreases  as  the  scatterer  radius  ka  is 
increased  from  0.5  (7.8%)  to  1.0  (7.0%). 

c.  For  both  cases,  the  maximum  error  occurs  at  kR=2.0. 

d.  If  the  maximum  error  points  (poles)  are  disregarded,  then 
the  maximum  error  is  less  than  or  equal  to  4%. 
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2.  For  the  dipolar  field: 

a.  The  normalized  error  at  kr  =  1.0  increases  by  a  factor  of  5, 
from  about  0.3%  to  about  1.5%  as  the  scatterer  radius  ka  is  increased  from  0.5 
to  1.0. 

b.  The  normalized  error  at  all  other  radii  is  essentially 
unchanged  by  varying  ka. 

c.  The  maximum  error  still  occurs  at  kr=2.0  as  it  occurred  for 
the  monopolar  field. 

d.  When  the  maximum  error  points  (poles)  are  disregarded 
then  the  maximum  error  is  less  than  3%. 

3.  For  the  quadrupolar  field: 

a.  The  normalized  error  at  a  given  node  shows  very  little 
dependence  on  the  value  of  ka. 

b.  The  maximum  error  remains  the  same  (4.1%)  as  the 
scatterer  radius  ka  is  increased  from  0.5  to  1 .0,  from  the  acoustic  center 

c.  For  both  cases,  if  the  maximum  error  points  are  disregarded 
then  the  error  is  less  than  2%. 

For  the  above  dipolar  and  quadrupolar  fields,  the  minimum  error  appears 
on  the  equatorial  points.  Moreover,  from  Figure  10,  very  similar  deviation  curves 
for  the  monopolar  (a,b),  dipolar  (c,d)  and  quadrupolar  (e,f)  fields,  for  ka=0.5  and 
ka=1.0  are  observed.  Also,  we  observe  that  the  maximum  normalized  error  on 
the  scatterer’s  surface  increases  as  the  scatterer’s  radius  increases. 

E.  INFLUENCE  OF  FLUID  MESH  ELEMENT  SIZE 

In  order  to  evaluate  the  influence  of  fluid  mesh  element  size  on  the 
results,  we  used  models  2  and  4.  Recall  that  model  2  is  divided  into  four  equal 
fluid  layers  of  thickness  L=0.5m  while  model  4  is  divided  into  two  equal  fluid 
layers  of  thickness  L=1.0m.  For  both  cases,  the  radiation  boundary  and 
scatterer’s  radius  satisfy  kR=2.5  and  ka=0.5,  with  k  =  1.0m'^ ,  R=2.5m  and 
a=0.5m. 
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Figure  11  presents  the  influence  of  the  element’s  size  (kL)  on  the 
normalized  scattered  pressure  deviation,  versus  kr.  The  monopolar  (a,b),  dipolar 
(c,d),  and  quadrupolar  (e,f)  fields  are  shown  from  top  to  bottom.  On  the  left  side 
the  element’s  size  is  kL=0.5  and  on  the  right  side  the  element’s  size  is  kL=1.0. 

Table  5  summarizes  the  results  for  the  maximum  error  in  percent  on  the 
scatterer’s  surface,  in  the  middle  and  at  the  end  of  each  layer  and  the  node 
location  for  the  monopolar,  dipolar  and  quadrupolar  fields,  when  the  element’s 
size  varies  and  kR,  ka  do  not. 


Type  of 

Field 

Scatterer 

Surface 

ka=0.5 

kr=1.0 

kr=1.5 

kr=2.0 

kr=2.5 

Maximum 

Normalized 

Error 

kr=1.5 

Monopolar 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

n=m=0 

kL=1.0 

1.2 

4.1 

6.3 

6.3 

4.5 

6.3  kr=2.0 

6.3 

kL=0.5 

1.7 

3.8 

5.7 

7.9 

3.8 

7.9  kr=2.0 

5.7 

Dipolar 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

n=1,  m=0 

kL=1.0 

0.7 

1.0 

2.6 

3.3 

3.9 

3.9  kr=2.5 

2.6 

kL=0.5 

0.6 

0.5 

2.2 

4.0 

3.6 

4.0  kr=2.0 

2.2 

Quadrupolar 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

n=2,  m=0 

kL=1.0 

0.2 

0.6 

0.9 

1.7 

4.2 

4.2  kr=2.5 

0.9 

kL=0.5 

0.2 

0.6 

0.7 

1.9 

4.0 

4.0  kr=2.5 

0.7 

Table  5.  Percent  normalized  scattered  pressure  deviation  and  location  of  the 
maximum  error,  for  monopolar,  dipolar  and  quadrupolar  incident  pressure  field, 
when  kL=1.0  and  kL=0.5,  scatterer  radius  ka=0.5,  radiation  boundary  kR=2.5. 


28 


Figure  11.  Normalized  scattered  pressure  deviation  computed  by  ATILA  versus 


kr.  From  top  to  bottom  monopolar  (a,b),  dipolar  (c,d)  and  quadrupolar  (e, 


incident  field;  on  the  left:  element  size  kL  =  0.5;  on  the  right;  element  size  kL 


1.0.  Scatterer  radius  a  and  radiation  boundary  R  satisfy:  ka  =  0.5  and  kR  =  2.5 


for  all  cases. 


II  O 


From  this  table  and  Figure  1 1  we  observe  that; 

1 .  For  the  monopolar  field: 

a.  On  the  scatterer’s  surface  ka=0.5,  as  the  element’s  size  is 
increased,  the  maximum  error  decreases  slightly  from  1.7%  to  1.2%. 

b.  The  maximum  error  at  each  layer  becomes  greater  as  kL  is 
increased,  except  at  the  maximum  deviation  point,  at  kr=2.0. 

c.  For  both  cases  the  maximum  error  occurs  at  kr=2.0. 

d.  Once  the  maximum  error  points  (poles)  are  disregarded, 
then  the  maximum  error  is  less  than  or  equal  to  4%.  Also,  in  that  case,  the 
corresponding  deviation  for  both  cases  is  similar. 

2.  For  the  dipolar  field: 

a.  On  the  scatterer’s  surface  for  both  cases  the  deviation  is 
almost  the  same. 

b.  As  the  element’s  size  is  increased,  the  error  increases 
except  at  the  maximum  deviation  point  at  kr=2.0. 

c.  The  maximum  deviation  for  kL=0.5  occurs  at  kr=2.0,  while 
for  kL=1.0  it  occurs  at  kr=2.5. 

d.  When  the  maximum  error  points  (poles)  are  disregarded, 
then  the  maximum  deviation  is  very  similar  and  less  than  3%  for  both  cases. 

3.  For  the  quadrupolar  field: 

a.  On  the  scatterer’s  surface,  ka=0.5,  the  maximum  deviation 
for  both  cases  is  0.2%. 

b.  For  both  cases,  the  variation  in  the  normalized  deviation 
follows  the  same  trend,  and  corresponding  values  are  very  close  to  each  other. 

c.  For  both  cases,  when  the  maximum  error  points  are 
disregarded,  the  deviation  is  less  than  2%. 

d.  The  maximum  deviation  occurs  at  kr=2.5  for  both  cases. 

F.  INFLUENCE  OF  FLUID  MESH  OUTER  RADIUS 

In  order  to  evaluate  the  influence  of  the  mesh  outer  radius,  we  used 
models  2  and  5.  Model  2  is  the  reference  model  and  is  divided  into  four  equal 
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fluid  layers  of  thickness  L=0.5m,  while  model  5  is  divided  into  two  equal  fluid 
layers  of  thickness  L=0.5m.  For  both  cases,  the  scatterer  radius  a  and  elements 
size  L,  satisfy  ka=0.5  and  kL=0.5,  for  k  =  1.0m"^  and  a=0.5m. 

Figure  12,  presents  the  influence  of  the  radiation  boundary  (kR)  on  the 
normalized  scattered  pressure  deviation  versus  kr.  The  monopolar  (a,b),  dipolar 
(c,d)  and  quadrupolar  (e,f)  fields  are  shown  from  top  to  bottom.  On  the  left  side 
the  radiation  boundary  is  kR=1.5  and  on  the  right  side  the  radiation  boundary  is 
kR=2.5. 

Table  6  summarizes  the  results  for  the  maximum  error  in  percent  in  the 
middle  of  each  layer  and  the  node  location  for  the  monopolar,  dipolar  and 
quadrupolar  fields  when  the  radiation  boundary  kR  varies  and  kL  and  ka  are  the 
same. 

It  is  observed  from  this  table  and  Figure  12  that : 

1 .  For  the  monopolar  field: 

a.  At  a  given  radius,  the  error  increases  as  the  radiation 
boundary  kR  is  increased.  For  example,  in  the  case  where  kR=2.5,  the  error  is 
100%  greater  than  the  corresponding  error  when  kR=1.5  at  the  location  where 
kr=1.5. 

b.  At  the  outer  fluid  boundary,  the  error  remains  almost  the 
same  (2%)  when  the  maximum  error  points  (poles)  are  disregarded. 

c.  When  kR=1.5,  the  maximum  error  occurs  at  the  outer  radius 
surface  points.  On  the  other  hand,  when  kR=2.5,  the  maximum  error  occurs  at 
kr=2.0. 

2.  For  the  dipolar  field: 

a.  The  error  decreases  drastically  as  we  increase  the  radiation 
boundary  kR  and  becomes  almost  100%  lower  at  the  boundary  for  kR=2.5. 

b.  The  maximum  error  still  occurs  at  the  same  radius  as  it  does 
for  the  monopolar  field. 
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2-EQUAL  FLUID  LAYERS  0.5M 


4-EQUAL  FLUID  LAYERS  0.5M 


Type  of 

Field 

Scatterer 

Surface 

ka=0.5 

Center  of 

Layer  a 

kr=0.75 

Center  of 

Layer  b 

kr=1.25 

Center  of 

Layer  c 

kr=1.75 

Center  of 

Layer  d 

kr=2.25 

Maximum 

Normalized 

Error 

kr=1.5 

Monopolar 

n=m=0 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

kR=1.5 

0.8 

1.4 

2.3 

2.3  kr=1.5 

2.3 

kR=2.5 

1.7 

2.9 

4.9 

6.9 

5.9 

7.8  kr=2.0 

5.7 

Dipolar 

n=1,  m=0 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

kR=1.5 

0.8 

1.6 

5.3 

— 

.... 

8.1  kr=1.5 

8.1 

kR=2.5 

0.5 

0.6 

1.8 

3.0 

3.9 

4.0  kr=2.0 

2.2 

Quadrupolar 

n=2,  m=0 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

Poles 

kR=1.5 

0.2 

0.2 

1.6 

.... 

.... 

4.8  kr=1.5 

4.8 

kR=2.5 

0.2 

0.3 

0.5 

1.2 

2.6 

4.1  kr=2.5 

0.7 

Table  6.  Percent  normalized  scattered  pressure  deviation  and  location  of  the 
maximum  error  for  monopolar,  dipolar  and  quadrupolar  incident  fields,  when 
kR=1 .5  and  kR=2.5,  scatterer  radius  ka=0.5  and  element  size  kL=0.5. 

c.  Once  again,  if  we  disregard  the  maximum  error  points 
(poles),  for  kR=2.5,  the  error  is  always  less  or  equal  to  3%.  On  the  other  hand, 
for  kR=1 .5,  the  error  is  less  than  7%. 

d.  The  minimum  error  points  for  both  cases  appear  on  the 

equator. 

3.  For  the  quadrupolar  field: 

a.  At  a  given  radius,  the  error  decreases  as  the  radiation 
boundary  increases. 

b.  The  maximum  error  for  kR=1.5  and  kR=2.5  occurs  on  the 
outer  boundary  surface. 

c.  For  both  cases,  if  the  maximum  error  points  (poles)  are 
disregarded  the  error  is  always  less  than  2%. 
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d.  Finally,  the  minimum  error  for  both  cases  appears  on  the 
equatorial  points. 

G.  SUMMARY  OF  RESULTS 

The  detailed  analyses  developed  for  the  five  models  listed  in  Chapter  III, 
Table  1,  indicate  that: 

1.  For  the  original  three  dimensional  model  (model  1),  the  maximum 
deviation  occurs  close  to  or  at  the  poles  for  the  axially  symmetric  monopolar 
(8.6%),  dipolar  (11.8%)  and  quadrupolar  (6%)  incident  pressure  fields.  Also,  the 
minimum  deviation  points  are  at  equatorial  nodes.  Moreover,  for  the  non-axially 
symmetric  incident  fields,  the  maximum  deviation  is  less  than  5.5%  and  occurs  at 
equatorial  points,  while  the  minimum  deviation  occurs  at  points  located  at  the 
poles. 

2.  Further  investigation  of  the  influence  of  the  fluid  mesh  inner  radius 
(ka),  the  element’s  size  (kL),  and  the  fluid  mesh  outer  radius  (kR),  on  the  results 
for  axially  symmetric  incident  pressure  fields,  revealed  that: 

a.  For  monopolar  and  dipolar  incident  fields,  when  the  poles 
are  disregarded: 

(1)  As  the  radiation  boundary  kR  increases  while  ka  and 
kL  remain  constant,  the  maximum  deviation  decreases  by  approximately  the 
same  amount  in  percent.  This  occurs  specifically  from  6.6%  when  kR=1.5,  to 
4%  when  kR=2.5,  for  ka=0.5  and  kL=0.5,  respectively,  which  corresponds  to 
65%. 

(2)  As  the  scatterer’s  radius  ka  increases,  while  kL  and 
kR  remain  constant,  the  maximum  deviation  decreases  by  a  small  amount  in 
percent.  This  occurs,  specifically,  from  4%  when  ka=0.5,  to  3.7%  when  ka=1.0, 
for  kL=0.5  and  kR=2.5,  respectively,  which  corresponds  to  8%. 

(3)  As  the  element’s  size  kL  increases,  while  ka  and  kR 
remain  constant,  the  maximum  deviation  remains  constant. 

b.  For  a  quadrupolar  incident  field,  the  maximum  deviation 
remains  essentially  the  same  between  1.9%  and  1.8%. 
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Table  7,  which  combines  three  tables,  summarizes  our  results  for  the 
maximum  normalized  error  in  percent,  when  the  poles  are  disregarded,  and  the 
type  of  field  for  which  this  error  occurs.  Also,  from  this  table  it  can  be  observed 
that  the  maximum  deviation  occurs  for  monopolar  and  dipolar  axially  symmetric 
incident  fields. 


kR=2.5 

kL=0.5 

kL=1.0 

ka=0.5 

4.0 

Monopolar 

4.0 

Dipolar 

ka=1.0 

3.7 

Monopolar 

ka=0.5 

kL=0.5 

kL=1.0 

kR 

6.6 

1.5 

Dipolar 

kR 

4.0 

4.0 

2.5 

Monopolar 

Dipolar 

kL=0.5 

ka=0.5 

ka=1.0 

kR 

6.6 

■KBjjjl 

1.5 

Dipolar 

kR 

4.0 

3.7 

2.5 

Monopolar 

Monopolar 

Table  7.  Percent  normalized  maximum  scattered  pressure  deviation  and  type  of 
incident  pressure  field,  when  poles  are  disregarded,  computed  by  ATILA. 
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V.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER 

INVESTIGATION 


A.  CONCLUSIONS 

Five  three-dimensional  spherical  models  were  developed  in  order  to 
evaluate  the  radiation  boundary  elements  used  in  the  ATILA  finite  element  code. 
The  models  simulate  a  rigid  spherical  solid  structure  surrounded  by  an  infinite 
fluid. 

Monopolar,  dipolar,  and  quadrupolar  incident  spherical  pressure  wave 
fields  were  imposed,  and  the  scattered  pressure  was  calculated  using  ATILA. 
Since  the  scattered  pressure  is  proportional  to  the  incident  pressure  and  has 
angular  dependence,  the  error  in  the  ATILA  results  at  each  finite-element  node 
was  quantified  by  normalizing  the  deviation  from  the  exact  value  by  the 
maximum  magnitude  of  the  scattered  pressure  at  the  same  radius.  This 
represents  the  normalized  scattered  pressure  deviation  computed  by  ATILA. 

The  range  of  values  of  the  dimensionless  lengths  which  characterized  the 
problem  was; 

The  fluid  mesh  inner  radius  ka  (scatterer  radius):  2.0 ,1.0, 0.5 
The  fluid  mesh  element  size  kL:  0.5 , 1 .0 
The  radiation  boundary  kR:  4.0 , 2.5 , 1 .5 

The  maximum  deviations  observed  were  for  the  axially  symmetric 
monopolar  (9%),  dipolar  (12%),  and  quadrupolar  (6%)  incident  pressure  fields. 
For  these  cases,  it  can  be  concluded  that  the  maximum  deviation  points  are  at 
the  poles,  while  the  minimum  deviation  occurs  on  the  equatorial  points. 

When  the  poles  are  disregarded,  for  the  monopolar  and  dipolar  incident 
pressure  fields  there  is  a  strong  influence  of  the  fluid  mesh  outer  radius  (kR)  and 
a  weak  influence  of  the  fluid  mesh  inner  radius  (ka)  on  the  results.  Specifically, 
as  the  fluid  mesh  outer  radius  increases  and/or  the  fluid  mesh  inner  radius 
increases  then  the  normalized  scattered  pressure  deviation  decreases.  For  the 
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quadrupolar  field,  the  maximum  deviation  remains  essentially  constant, 
independent  of  the  above  factors  (ka,  kR).  Furthermore,  variation  of  the 
elements  size  (kL)  does  not  affect  the  maximum  normalized  scattered  pressure 
deviation. 

Moreover,  for  the  non-axially  symmetric  incident  pressure  fields,  the 
maximum  deviation  observed  was  less  than  5.5%,  and  was  found  to  occur  at 
nodes  on  the  equator,  while  the  minimum  deviation  points  were  located  at  the 
poles. 

B.  SUGGESTIONS  FOR  FURTHER  INVESTIGATION 

The  magnitudes  of  the  errors  in  the  scattered  pressures  computed  by 
ATI  LA  were  larger  than  expected,  based  upon  the  mesh  sizes  employed  in  the 
models;  for  all  of  the  models  the  radial  length  L  of  the  fluid  elements  was  less 
than  1/6  wavelength.  It  would  have  been  desirable  to  examine  the  calculated 
scattering  for  a  more  refined  mesh.  However,  to  divide  the  element  scale  size  by 
a  factor  of  2  would  have  resulted  in  models  which  exceed  the  allowed  number  of 
degrees  of  freedom.  Of  course,  this  is  because  all  the  models  employed  were 
fully  three-dimensional,  and  did  not  make  use  of  any  possible  reductions  in  size 
due  to  wave  function  symmetry  (in  fact,  our  choice  of  e'*  for  the  azimuthal  wave 
function  precluded  reduction  in  that  direction).  Also,  we  are  interested  in  the 
computed  scattering  in  all  of  the  spherical  harmonic  wave  components  for  an 
incident  wave  of  only  one  component. 

It  is  suggested  that  further  investigation  of  the  performance  of  the  ATILA 
radiation  boundary  elements  be  conducted  using  two-dimensional  models.  It 
might  also  be  advantageous  to  use  sin  ^  or  cos  <j)  instead  of  to  represent  the 
azimuthal  dependence.  This  will  allow  a  finer  mesh  to  be  employed,  and  a 
broader  range  of  values  of  ka,  kL,  and  kR  to  be  examined. 
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APPENDIX  A.  C  CODE  FOR  ANALYTICAL  PRESSURE 
SCATTERED  CALCULATION 


************************************************************************ 

*  PROGRAM  SCHAR.C  ,  by  Panagiotis  Sinanoglou  2/2/1996 

*  Program:Computes  the  Scattered  Spherical  Harmonic  pressure  from  a  rigid 

*  spherical  structure  at  the  point  r1  (x,y,z),  for  the  Wavenumber  k,  for; 

*  natural  boundary  and  pressure  release  condition 

*  Input  Variables: 

*  x,y,z;  Cartesian  Coordinates  of  the  Node  Point. 

*  k:  Wavenumber. 

*  r2:  Radius  of  the  Scatterer's  Surface. 

n,m:  Orders  of  Spherical  Hankel  and  Legendre  functions. 

*  Output  Variables: 

*  r1 :  Radius  in  meters. 

*  phi:  Azimouthial  angle  in  degrees. 

*  theta:  Polar  angle  in  degrees. 

*  rp5:  Real  part  of  computed  scattered  pressure  at  r1  in  pascals 

*  rp6:  Imaginary  part  of  computed  scattered  pressure  at  r1  in  pascals 

*  Formula  for  Incident  Wave:  Pinc=Pnm(costheta)*e''(imphi)*hn(1)(kr2) 

*  Formula  for  Scattered  Wave  from  Rigid  Boundary: 

*  Psc=-Pnm(costheta)*e''(imphi)*[hn'(1  )(kr2)/hn'(2)(kr2)rhn2(kr1 ) 

*  Formula  for  Scattered  Wave  from  Pressure  Release  Surface: 

*  Psc=-Pnm(costheta)*e'^(imphi)*[hn(1  )(kr2)/hn(2)(kr2)]*hn2(kr1 ) 

*  External  Functions:  sphbes.c,  pigndr.c,  bessjy.c,  beschb.c,  chebev.c, 
complex.c,  nrutil.c 

*  Header  Files  for  prototype  function  declaration:nr.h  nrutil.h.complexl.h 

*  Functions: 

*  Normalized  Legendre[  plgndr(n,m,x)  ], Spherical  Hankel  for  the 

*  real  and  the  imaginary  parts[  sphbes(n,(k*r2),&xsj.&xsy,&xsjp,&xsyp)  ], 

*  Complex(structures)[  Cexp(mphi),  Cdiv(a,b),  Cmul(a,b),  RCmul(a,b)  ] 

*  File  nphraSSn.c  provides  node  coordinates(x,y,z). 

*  File  diml.dt  stores  the  corresponding  spherical  coordinates(r1, theta, phi) 

*  File  dimIS.dt  stores  the  Legendre  plgndr(n,m,(z/r1). 

*  File  dimi  1  .dt  stores  the  real, imaginary  part  of  exp.  function  e''(imphi). 

*  File  dim14.dt  stores  the  ratio  of  the  derivatives  of  the  Spherical  Hankel  for  the 
real  and 

*  imaginary  part  of  hn'(1)(kr2)/hn'(2)(kr2)  for  natural  boundary  conditions  and  the 
ratio  of 

*hn(1)(kr2)/hn(2)(kr2)  for  pressure  release,evaluated  at  the  scattering  surface,  at 
r2=0.5m. 

*  File  dimIS.dt  stores  the  real  and  the  imaginary  part  of  the  Spherical 

*  Hankel  hn2(kr1 )  at  the  node  range  r1 . 
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*  File  dim16.dt  stores  the  real  and  the  imaginary  part  of  the  scattered  pressure 
Psc,  at  the  *  range  r1 . 

*********************************************************************** 


#include  <stdio.h> 

#include  <math.h> 

#include<string.h> 

#include<stdlib.h> 

#include  "complexi  .h" 

#define  NRANSI 
#include  "nr.h" 

#include  "nrutil.h" 

#define  pi  3.141592654 

main() 

{ 

FILE  *f1 ;  FILE  *f3;  FILE  *f2:  FILE  *f4;  FILE  *f5;FILE  *f6;FILE  ‘fT; 
FILE  *f20:  FILE  *f21;FILE  *f22; 
fcomplex  a,b,hn2; 
float  r2=0.5,k=2.0; 

float  sj,sy,sjp,syp,xsj,xsy,xsjp,xsyp,rp1  ,rp2,rp3,rp4,rp5,rp6; 
intj.rpm; 

float  fac,val,t,f1 1  ,f12,f13,f14,Pnm,s1  ,s2; 
unsigned  int  factorial(unsigned  int  a); 
double  x,y,z,z1  ,r,r1  .theta, phi, q.mphi; 
int  d,i,N,ch,m,n; 

charg[11].h[11],c[11],e[11].f[11]: 

f3=fopen("dim1.dt","w"); 
printffEnter  the  integer  valuesfor  N,n,m\n"); 
scanf("%d  %d  %d",&N,&n,&m): 
f4=fopen("sphra33b.c","r"): 
f5=fopenrdim11.dr,"w"); 
f6=fopenrdim12.dt"."w"): 
f7=fopen("dim13.dt",'W'): 
f20=fopen("dim14.dt","w"): 
f21=fopen("dim15.dt","w"): 
f22=fopen("dim16.dt","w"): 
for(i=0:  i<N;  ++i) 

{ 

fscanf(f4,"%lf  %lf  %lf  %d 

%lf%s%s%s%s%s\n",&x,&y,&z,&d,&q,&g,&h,&c,&e,&f); 

/*  Converts  from  cartesian(x,y,z)  to  spherical  coordinates(r1,theta,phi).7 
if((x>0.0  &&  y>=0.0)) 

{ r1=sqrt(x*x+y*y+z*z): 
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phi=atan(y/x); 

theta=acos(z/r1);} 

else  if(x==0.0  &&  y==0.0  &&  z==0.0) 

{r1=0.0;theta=0.0:  phi=0.0;} 
else  if  (x==0.0  &&  y==0.0  &&  z!=0.0) 

{ r1=sqrt(x*x+y*y+z*z); 
phi=0.0: 

r  1  =sq  rt(x*x+y*y+z*z) ; 
theta=acos(z/r1);} 

else  if(x==0.0  &&  y>0.0  &&(z==0.0|lz!=0.0)) 

{ phi=pi/2: 

r1  =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1):} 
else  if(y==0.0  &&  x<0.0  ) 

{  phi=pi: 

r1  =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 
else  if(y<0.0  &&  x==0.0) 

{  phi=3*(pi/2); 
r1  =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 
else  if  ((x<0.0  &&  y>0.0) ) 

{phi=+(pi)+atan(y/x): 
r1  =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1); } 
else  if  (x<0.0  &&  y<0.0) 

{  phi=pi+atan(y/x): 
r1  =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 
else  if  (x>0.0  &&  y<0.0) 

{  phi=(2*pi)+atan(y/x): 
r1  =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1):} 

theta=theta*(1 80/pi); 
phi=phi*360/(2*pi): 

fprintf(f3,"  %+4.3lf  %+6.3lf  %+6.3lf\n",r1. (theta), (phi)); 

/*  The  required  value  of  xx=cos(theta)  for  the  Pnm(xx)  is:  (z/r1)  */ 

/*  Calculates  the  Normalized  Legendre  function  "plgndr(n,m,(z/r1)  */ 
fprintf(f7,"%4s  %4s  %10s  %24s\n","n","m","x","plgndr(n,m,x)"); 
rpm=abs(m); 

if(m>=0  ) 

{ f1 3=factorial(n-m); 
f14=factorial(n+m); 
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s1=sqrt((2*n+1)*f13/(4*pm4)); 

Pnm=(plgndr(n,m,(z/r1))*s1): 

fprintf(f7,"%4d  %4d  %13.6lf  %19.6e\n".n,m,(z/r1),Pnm): 

} 

if(m<0) 

{ f1 3=factorial(n-rpm): 
f14=factorial(n+rpm); 
s2=sqrt((2*n+1)*f14/(4*pi*f1 3)); 
f1 1=factorial(n-rpm); 
f1 2=factorial(n+rpm); 
r=pow(-1  ,rpm): 
t=r*(f11/f12); 

Pnm=(plgndr(n,rpm,(z/r1))*t*s2); 
fprintf(f7,"%4d  %4d  %13.6lf  %19.6e\n",n,m,(z/r1),Pnm): 

} 

/*  Calculates  the  e''(imphi)  Real/lmagin  parts  */ 
mphi=(m*phi)*2*{pi/360); 

fprintf(f6,"%lf  %lf  \n",Cexp(mphi).r,Cexp(mphi).i): 

/*  Calculates  the  (z/r1),phi(degrees)  */ 

fprintf(f5,"%+6.3lf  %+6.3lf  \n",(z/r1),(phi)): 

/*  Calculates  the  ratio  of  the  Spherical  Hankel  for  the  real  and  the  imaginary 
parts  */ 

/*  of  hn'(1)(kr2)/hn'(2)(kr2)  on  the  scattering  surface  at  r2=0.5  for  N.B.C  ,  or  the 
*/ 

/*  ratio  of  hn(1)(kr2)/hn(2)(kr2)  for  pressure  release  condition 
*/ 


sphbes(n,(k*r2),&xsj,&xsy,&xsjp,&xsyp): 
a.r=xsjp;  /*  a.r=xsj;  */ 

a. i=xsyp;  /*  a.i=xsy:  */ 

b. r=xsjp;  /*  b.r=xsj;  */ 

b.i=-xsyp:  /*  b.i=-xsy;  */ 

fprintf(f20,"%30s  \n","The  real  and  the  imag.  of  hn'(1)/hn'(2)  is:"); 
fprintf(f20,"%f  %f  \n",Cdiv(a,b).r,Cdiv(a,b).i); 

/*  Calculates  The  real  and  the  imag.  of  hn2(kr1)  */ 
sphbes(n,(k*r1),&xsj,&xsy,&xsjp,&xsyp); 
hn2.r=xsj: 
hn2.i=-xsy: 

fprintf(f21,"%30s  \n","The  real  and  the  imag.  of  hn2(r1)  is:"); 
fprintf(f21  ,"%f  %f  \n",hn2.r,hn2.i): 

/*  Calculates  The  real  and  the  imag.parts  of  the  scattered  pressure  at  the  */ 

/*  radius  (r1)  :Psc=-Pnm(costheta)*e'^(imphi)*[hn'(1)(kr2)/hn'(2)(kr2)]*hn2(kr1)*/ 
rp5=Cmul(RCmul(-Pnm,Cmul(Cdiv(a,b),Cexp(mphi))),hn2).r; 
rp6=Cmul(RCmul(-Pnm,Cmul(Cdiv(a,b),Cexp(mphi))),hn2).i; 
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fprintf(f22,"%f  %f\n",rp5,rp6); 

} 

fclose(f4): 

fclose(f3); 

fclose(f5): 

fclose(f6): 

fclose(f7); 

fclose(f20); 

fclose(f21): 

fclose(f22): 

return  0; 

} 

#undef  NRANSI 


unsigned  int  factorial(unsigned  int  a) 

{ 

if((a==1)||(a==0)) 
return  1; 
else 

{a*=factorial(a-1); 
return  a; 

} 

} 

************************************************************************ 
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APPENDIX  B.  FUNCTION  INCPRE(X,Y,Z,K) 


FUNCTION  INCPRE(X,Y,Z,K) 


*  PROGRAM  BY  ARTHUR  LOBO  DA  COSTA  RUIZ  12/23/93 

*  MODIFIED  BY  STEVEN  R.  BAKER  3/4/96 


*  FUNCTION: 

*  COMPUTES  THE  INCIDENT  SPHERICAL  HARMONIC  PRESSURE  AT  THE 
‘POINT  (X.Y,Z),FOR  THE  WAVENUMBER  K 

*  VARIABLES  INPUT: 

*  X,Y.Z:  CARTESIAN  COORDINATES  OF  THE  POINT 

*  K:  WAVENUMBER 

*  VARIABLES  OUTPUT: 

*  RADIUS:  R  IN  METERS 

*  PHI  IN  DEGREES  (HORIZONTAL  ANGLE) 

*  THETA  IN  RADIANS  (AZIMUTAL  ANGLE) 

*  INCPRE  :  PRESSURE  AT  POINT  IN  PASCAL 

* 

*  FORMULA  USED: 

* 

*  INCPPRE=  H1(N)*P(M.N)*DCMPLX(DCOSD(M*PHI),DSIND(M*PHI)) 

*  THIS  MAKES  AN  OUTGOING  INCIDENT  WAVE  FOR  E'^IWT  TIME 

*  DEPENDENCE  ,  AS  FOR  ATILA 


DOUBLE  PRECISION  K.X.Y.Z.R, PHI, THETA, KR 
REAL*8  F1,PMN(-2:2,0:2),LOUT 
INTEGER  N,M,NMAX 
COMPLEX*16  INCPRE,H1(0:2),H1OUT 

*  N  AND  M  ARE  ORDERS  FOR  HANKEL  AND  LEGENDRE  FUNCTIONS 

N=0 

M=0 

*  HI  OUT  AND  LOUT  ARE  HANKEL  AND  LEGENDRE  OUTPUTS 

*  REF  TON  AND  M 

NMAX=2 

*  TRANSFORM  CARTESIAN  COORDINATES  (X,Y,Z) 

*  INTO  SPHERICAL  COORDINATES 

*  R(RADIUS),  PHI(/\ZIMUTAL  ANGLE)  AND  THETA(POLAR  ANGLE) 

R=DSQRT(X*X+Y*Y+Z*Z) 

PHI=DATAN2D(Y,X) 

IF((X.EQ.O).AND.(Y.EQ.O))  PHNO.ODO 
THETA=DACOS(Z/R) 
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KR=K*R 

NMAX  IS  THE  MAXIMUM  NUMBER  OF  HARMONICS 
CALL  HANKEL1(KR,NMAX,H1) 

SUBROUTINE  HANKEL1  RETURNS  SPHERICAL  HANKEL 

FUNCTIONS  OF  THE  FIRST  KIND  JN  +  I  YN 
CALL  LEGNDR(THETA,NMAX,PMN) 

LOUT=F1(N,M)*PMN(M,N) 

SUBROUTINE  LEGNDR  RETURNS  ASSOCIATE  LEGENDRE  FUNTION 
H10UT=H1(N) 

INCPRE=H10UT*L0UT*DCMPLX(DC0SD(M*PHI),DSIND(M*PHI)) 

IF  ((R.LE.0.501).AND.(R.GT.0.499))  THEN 
PRINT  *,"X,Y,Z  =",X,Y,Z 
PRINT  *, PHI, THETA 
PRINT  V'K.R.KR  =".K,R,KR 
PRINT  *,"REAL  HANKEL1  =",H10UT 
PRINT  *, "LEGENDRE  =",LOUT 
PRINT  MNCPRE 
PRINT  *  ”***************" 

ELSE 

CONTINUE 

ENDIF 

RETURN 

END 

‘  I**************************************************************** 

SUBROUTINE  HANKEL1(X,NMAX,H1) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEXES  H1(0:NMAX) 


C  GIVEN  THE  VARIABLE  X,  AND  THE  MAXIMUM  ORDER  NMAX, 

C  THIS  ROUTINE  GENERATES  THE  SPHERICAL  HANKEL  FUNCTION  OF 
THE 

C  FIRST  KIND  H1N  FOR  ALL  N  FROM  0  TO  NMAX  (INCLUSIVE) 

C  INPUT: 

C  X  =  DOUBLE  PREC.  VARIABLE  (RADIUS) 

C  NMAX  =  INTEGER  MAXIMUM  ORDER  OF  BESSEL  FUNCTIONS 
DESIRED 
C  OUTPUT: 

C  H1(N)  =  ARRAY  OF  SPHERICAL  HANKEL  FUNCTIONS  H1N(X),  WHERE 
C  H1N  =  JN  +  IYN 

C  THIS  ROUTINE  IS  BASED  ON  THE  RECURSION  FORMULA 
C  FROM  ABRAMOWITZ  &  STEGUN:  10.1.10  &  10.1.15,  PP.438-9 
C  THE  PS  ARE  THE  COEFFICIENTS  OF  ORDER  N  &  -(N+1 ), 

C  THE  FO'S  ARE  OLD  PS,  FOR  RECURSION 
IF  ( X  .LE.  O.ODO  )  THEN 
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H1(0)  =  DCMPLX(1.0D0,-1.0D35) 

D02  N  =  1,  NMAX 
H1(N)  =  CMPLX(0.0D0,-1.0D35) 

2  CONTINUE 
RETURN 
END  IF 
SX  =  DSIN(X) 

CX  =  DCOS(X) 

XINV=  1.0D0/X 
M1N  =  -1.0D0 
FN  =  XINV 
FMN  =  O.ODO 
FNO  =  FMN 
FMNO  =  FN 
DO  4  N  =  0,  NMAX 

H1(N)  =  CMPLX(  FN*SX  +  M1N*FMN*CX,  -FN*CX  +  M1N*FMN*SX  ) 

T1  =  (2*N+1)*XINV 
T2  =  T1*FN-FNO 
FNO  =  FN 
FN  =  T2 

T2  =  -T1*FMN-FMNO 
FMNO  =  FMN 
FMN  =  T2 
M1N  =  -M1N 
4  CONTINUE 
RETURN 
END 

0  *ic’:kie**icif-k***ie*i('kic*icieic***itif’ki(’kic*icic'k-k1fic-k*ieie*i(icic*-kic‘kie*ie'kiei:ieie-k-k*-k*-k* 

SUBROUTINE  LEGNDR(THETA,NMAX,PMN) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  F1,PMN(-NMAX:NMAX,0:NMAX) 

C  GIVEN  THE  VARIABLE  THETA,  AND  THE  MAXIMUM  ORDER  NMAX, 

C  THIS  ROUTINE  GENERATES  THE  ASSOC.  LEGENDRE  FUNCTIONS  PMN 
C  OF  THE  ARGUMENT  COS(THETA)  (THETA  MUST  BE  BETWEEN  0  &  PI) 
C  FOR  ALL  N  FROM  0  TO  NMAX  (INCLUSIVE) 

C  AND  FOR  ALL  M  FROM  -N  TO  N  (SOME  OTHERS  SET  TO  ZERO) 

C  INPUT: 

C  THETA  =  VARIABLE  (POLAR  ANGLE),  MUST  BE  BETWEEN  0  &  PI 
(INCL) 

C  NMAX  =  INTEGER  MAXIMUM  ORDER  OF  LEGENDRE  FUNCTIONS 
C  DESIRED 
C  OUTPUT: 

C  PMN  =  DOUBLE  PREC.  ARRAY,  CONTAINS  ASSOC.  LEGENDRE  FNS 
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C  THIS  ROUTINE  IS  BASED  ON  THE  RECURSION  FORMULAE 

C  FROM  ABRAMOWITZ  &  STEGUN 
X  =  DCOS(THETA) 

SINTHT  =  DSIN(THETA) 

IF  (  SINTHT  .GT.  0.  )  THEN 
S1NINV=  1.0D0/SINTHT 
ELSE 

SININV  =  0.0D0 
END  IF 

C  SET  VALUES  FOR  N  =  0,  1  (NMAX  MUST  BE  AT  LEAST  1) 

PMN{0,0)  =  1.0D0 
PMN{1,0)  =  0.0D0 
PMN(-1,0)  =  0.0D0 
PMN(0,1)  =  X 
PMN{1,1)  = -SINTHT 
PMN(-1,1)  =  SINTHT*0.5D0 

C  IN  LOOP,  TNP1  =  2*N+1,  TNP2FC  =  (2*N+2)!,  M1N  =  (-1)**(N+1) 
TNP1  =  1.0D0 
TNP2FC  =  2.0D0 
M1N  =  -1.0D0 
DO  4  N  =  1 ,  NMAX-1 
TNP1  =  TNP1  +  2.0D0 
TNP2FC  =  TNP2FC  *  TNP1  *  (TNP1+1) 

M1N=-M1N 
DO  3  M  =  -N.  N 

PMN(M,N+1)  =  (TNP1*X*PMN(M,N)  -  (N+M)*PMN(M,N-1))/(N-M+1) 

3  CONTINUE 
PMN(N+1.N)  =  0.0D0 
PMN(-N-1,N)  =  0.0D0 

PMN(N+1,N+1)  =  (X*PMN(N,N+1)  -TNP1*PMN(N,N))  *  SININV 
PMN(-N-1,N+1)  =  M1N*PMN(N+1,N+1)/TNP2FC 

4  CONTINUE 
RETURN 

END 

******************************* 

FUNCTION  F1(NN.MM) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  FACT1,FACT2,PI 
PI=4.0D0*DATAN(1.0D0) 

FACT1=1.0D0 
DO  10  l=1,NN+MM 
10  FACT1=FACT1*I 
FACT2=1.0D0 
DO  20  l=1,NN-MM 


48 


20  FACT2=FACT2*I 

F1=DSQRT((2*NN+1)*FACT2/4.0D0/PI/FACT1) 

RETURN 

END 
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